void calc_orb(P *p, TRHO *rho, int *J, int nn, bool *failflag ) {
	
			if (8<nn) { 
				printf("ERROR: nn=%i, max is 8\n", nn);
				exit(0);
			}
	
			REAL pos[3][8], vel[3][8], dt[8], dt_max[8];
			REAL acc[3][8], t[8], r[8], z[8], R[8], ar[8], az[8], a_mod[8], dt_new[8], _dt[8];
			
			INT norbit[8];
			REAL radsign[8], radsign_previous[8];
			FILE *__fp[8];

			
			bool resultF = true;
			int nsteps = ns*np;
			REAL ins2PI = 1.0/ns*2*M_PI;
			REAL eps2 = SQR(eps);
			
			__m128 ins2PIjj = _mm_set1_ps(ins2PI);
			
			for (int j=0; j<nn; j++)
			memset(&rho[J[j]*dg.n], 0, dg.n*sizeof(*rho));
			
			
			for (int j=0; j<nn; j++) {
				
			//memcpy(pos[j], p[J[j]].Pos, 3*sizeof(REAL)); 
			//memcpy(vel[j], p[J[j]].Vel, 3*sizeof(REAL));
			
			pos[0][j] = p[J[j]].Pos[0];
			pos[1][j] = p[J[j]].Pos[1];
			pos[2][j] = p[J[j]].Pos[2];
			
			vel[0][j] = p[J[j]].Vel[0];
			vel[1][j] = p[J[j]].Vel[1];
			vel[2][j] = p[J[j]].Vel[2];
			
			
			
			REAL r2 = SQR(pos[0][j])+SQR(pos[1][j]);
			REAL r = sqrt(r2);
			REAL z = fabs(pos[2][j]);
			
			// serial
			
			AG::A a = ag.get_acc_pot(r,z);
			
			
			
			REAL a_mod = sqrt( SQR(a.r)+SQR(a.z) );
			
			if (r!=0) {
				REAL a_r_over_r(a.r/r);
				acc[0][j]= pos[0][j]*a_r_over_r;
				acc[1][j]= pos[1][j]*a_r_over_r;
			} 
			
			if (z!=0)
				acc[2][j]= pos[2][j]*a.z/z;
	
			
			REAL posm = sqrt( SQR(pos[0][j]) + SQR(pos[1][j]) + SQR(pos[2][j]) );
			REAL velm = sqrt( SQR(vel[0][j]) + SQR(vel[1][j]) + SQR(vel[2][j]) );
			REAL accm = sqrt( SQR(acc[0][j]) + SQR(acc[1][j]) + SQR(acc[2][j]) );
			
			
			REAL tp = 2*M_PI*posm/sqrt(accm*posm);
			dt[j] = tp/ns;

			
			t[j] = 0;
			
			
			dt_max[j] = 0.1*2*M_PI*r*sqrt(-a.pot);
			
			// 0. predict velocity at dt/2
			vel[0][j] += acc[0][j]*.5*dt[j];
			vel[1][j] += acc[1][j]*.5*dt[j];
			vel[2][j] += acc[2][j]*.5*dt[j];
		
			failflag[j] = false;
			
			norbit[j]=0;
			radsign_previous[j] = radsign[j] = pos[0][j]*vel[0][j] + pos[1][j]*vel[1][j] + pos[2][j]*vel[2][j];
			
			//char fname[256];
			//sprintf(fname, "orbit.%i.txt", j);
			//__fp[j] =  fopen(fname, "w");

			
			}
			
			static const __m128 sign_mask = _mm_set1_ps(-0.f); // -0.f = 1 << 31
			static const __m128 nulljj = _mm_set1_ps(0);
			//return _mm_andnot_ps(sign_mask, x);

			
			__m128 pos0jj = _mm_load_ps(pos[0]);
			__m128 pos1jj = _mm_load_ps(pos[1]); 
			__m128 pos2jj = _mm_load_ps(pos[2]);
				
			__m128 vel0jj = _mm_load_ps(vel[0]);
			__m128 vel1jj = _mm_load_ps(vel[1]); 
			__m128 vel2jj = _mm_load_ps(vel[2]); 
			
			
			//__m128 failflagjj = _mm_set1_ps(0);
			
 			for (int i=0; i<100*nsteps; i++) {
				
				__m128 dtjj =_mm_load_ps(dt);
				
				
				pos0jj =_mm_add_ps(pos0jj, _mm_mul_ps(vel0jj, dtjj));
				pos1jj =_mm_add_ps(pos1jj, _mm_mul_ps(vel1jj, dtjj));
				pos2jj =_mm_add_ps(pos2jj, _mm_mul_ps(vel2jj, dtjj));
				
				__m128 r2jj = _mm_add_ps(_mm_mul_ps(pos0jj, pos0jj), _mm_mul_ps(pos1jj, pos1jj)); 
				__m128 rjj = _mm_sqrt_ps(r2jj);
				
				__m128 zjj = _mm_andnot_ps(sign_mask, pos2jj);
				
				__m128 R2jj = _mm_add_ps(r2jj, _mm_mul_ps(pos2jj, pos2jj)); 
				__m128 Rjj = _mm_sqrt_ps(R2jj);

				_mm_store_ps(r, rjj);
				_mm_store_ps(z, zjj);
				//_mm_store_ps(R, Rjj);

				_mm_store_ps(pos[0], pos0jj);
				_mm_store_ps(pos[1], pos1jj);
				_mm_store_ps(pos[2], pos2jj);

				
				__m128 arjj, azjj;
				ag.get_acc_ps(rjj,zjj, arjj, azjj);

				//failflagjj = _mm_or_ps(failflagjj, _mm_and_ps(_mm_cmpeq_ps(arjj, nulljj), _mm_cmpeq_ps(arjj, nulljj)));

				_mm_store_ps(ar, arjj);
				_mm_store_ps(az, azjj);
				

				
				////  REAL a_mod = sqrt(SQR(ar[jj])+SQR(az[jj]));  ////
				////  REAL dt_new = ins2PI*R[jj]/sqrt(a_mod*R[jj]); ////
				
				//REAL a_modR[4];
				__m128 a_modjj = _mm_sqrt_ps(_mm_add_ps(_mm_mul_ps(arjj,arjj), _mm_mul_ps(azjj,azjj)));
				__m128 dt_newjj = _mm_div_ps(_mm_mul_ps(ins2PIjj, Rjj), _mm_sqrt_ps(_mm_mul_ps(a_modjj, Rjj)));
				_mm_store_ps(dt_new, dt_newjj);
				
				
				
				
				
				
			
				/*
				// 1. predict new position 
				pos[0][jj] += vel[0][jj]*dt[jj];
				pos[1][jj] += vel[1][jj]*dt[jj];
				pos[2][jj] += vel[2][jj]*dt[jj];
				
				
				
				// 2. get accelerations
				REAL r2 = SQR(pos[0][jj])+SQR(pos[1][jj]);
				r[jj] = sqrt(r2);
				z[jj] = fabs(pos[2][jj]);
				
				REAL R2 = r2 + SQR(pos[2][jj]);
				R[jj] = sqrt(R2);
				*/
				/*
				AG::A  a = ag.get_acc(r[jj],z[jj]);
				
				ar[jj] = a.r;
				az[jj] = a.z;
				*/
					
				//REAL a_mod = sqrt(SQR(ar[jj])+SQR(az[jj]));
				//REAL dt_new = ins2PI*R[jj]/sqrt(a_mod[jj]*R[jj]);

				
				
				
				// 2. get accelerations
				if (dg_dtmax_cell)
				for (int jj=0; jj<nn; jj++) {
				#ifdef DG_DTMAX_CELL
				REAL vr = fabs(vel[0][jj]*pos[0][jj] + vel[1][jj]*pos[1][jj])/r[jj];
				REAL vz = fabs(vel[2][jj]);
				
				int ir = dg.r2ir(r[jj]);
				int iz = dg.z2iz(z[jj]);
				
				REAL dr = dg.ir2r(ir+1)-dg.ir2r(ir);
				REAL dz = dg.iz2z(iz+1)-dg.iz2z(iz);
				
				REAL dtrmax = dr/vr;
				REAL dtzmax = dz/vz;
				
				if (dtrmax<dt_new[jj]) dt_new[jj] = dtrmax; 
				if (dtzmax<dt_new[jj]) dt_new[jj] = dtzmax;
				#endif
				}
				
				
				
				int n_finished(0);
				for (int jj=0; jj<nn; jj++) {
				
					
					if (norbit[jj]<np) {

						if (ar[jj]==0 && az[jj]==0) {
							n_finished++;
							failflag[jj]=true;
						}
						
						radsign[jj] = pos[0][jj]*vel[0][jj] + pos[1][jj]*vel[1][jj] + pos[2][jj]*vel[2][jj];

						if(radsign[jj] > 0 && radsign_previous[jj] < 0) norbit[jj]++;
						radsign_previous[jj] = radsign[jj];

						
						if (dt_max[jj]<dt_new[jj]) dt_new[jj] = dt_max[jj];
				
						_dt[jj] = 0.5*(dt[jj] + dt_new[jj]);
			
						t[jj] +=_dt[jj];
						dt[jj] = dt_new[jj];
					
						#ifdef DG_ADD_STD
						dg.add(&rho[J[jj]*dg.n], r[jj], z[jj], _dt[jj]);
						#endif
					
						//fprintf(__fp[jj], "%g %g %g\n", pos[0][jj], pos[1][jj], pos[2][jj]);
					} else n_finished++;
				}	
				
				if (n_finished==nn) break;
				
				
				
				/*
				#ifdef DG_ADD_STD
				dg.add_ps(rho, J, rjj, zjj, _dt);
				#endif
				*/
				
				//dg.add(&rho[J[jj]*dg.n], r[jj], z[jj], _dt);
				
				/*
			}
			
				//if (!failflag[j]) 
			#ifdef DG_ADD_STD
			dg.add_ps(rho, J, rjj, zjj, _dt);
			#endif
			)
			for (int jj=0; jj<nn; jj++) {		
			*/
				
				
				/*
				if (r!=0) {
					REAL a_r_over_r(ar[jj]/r[jj]);
					acc[0][jj]= pos[0][jj]*a_r_over_r;
					acc[1][jj]= pos[1][jj]*a_r_over_r;
				} 
			
				if (z!=0)
 					acc[2][jj]= (pos[2][jj]*(az[jj])/z[jj]);
				
				
				vel[0][jj] += acc[0][jj]*_dt[jj];
				vel[1][jj] += acc[1][jj]*_dt[jj];
				vel[2][jj] += acc[2][jj]*_dt[jj];
				*/
				/*
				if (jj==0) {
					//a.print();
					printf("%4i % .7e, % .7e\n", i, r[0], z[0]);
				//if (r[0]!=r[0]) exit(0);
				}
					*/
			
			
			__m128 a_r_over_rjj = _mm_div_ps(arjj, rjj);
			__m128 a_z_over_zjj = _mm_div_ps(azjj, zjj);
			
			__m128 maskrjj = _mm_cmpneq_ps(rjj, nulljj);
			__m128 maskzjj = _mm_cmpneq_ps(zjj, nulljj);
			
			a_r_over_rjj = _mm_and_ps(maskrjj, a_r_over_rjj);
			a_z_over_zjj = _mm_and_ps(maskzjj, a_z_over_zjj);			
			
			__m128 acc0jj =_mm_mul_ps(pos0jj, a_r_over_rjj);
			__m128 acc1jj =_mm_mul_ps(pos1jj, a_r_over_rjj);
			__m128 acc2jj =_mm_mul_ps(pos2jj, a_z_over_zjj);
			/*
			_mm_store_ps(acc[0], acc0jj);
			_mm_store_ps(acc[1], acc1jj);
			_mm_store_ps(acc[2], acc2jj);
			/**/
			/*
			for (int jj=0; jj<nn; jj++) {
				// 3. predic new velocities				
				vel[0][jj] += acc[0][jj]*_dt[jj];
				vel[1][jj] += acc[1][jj]*_dt[jj];
				vel[2][jj] += acc[2][jj]*_dt[jj];
			}
			*/
			
			__m128 _dtjj = _mm_load_ps(_dt);
			vel0jj =_mm_add_ps(vel0jj, _mm_mul_ps(acc0jj, _dtjj));
			vel1jj =_mm_add_ps(vel1jj, _mm_mul_ps(acc1jj, _dtjj));
			vel2jj =_mm_add_ps(vel2jj, _mm_mul_ps(acc2jj, _dtjj));

			_mm_store_ps(vel[0], vel0jj);
			_mm_store_ps(vel[1], vel1jj);
			_mm_store_ps(vel[2], vel2jj);
			
			
			
				//printf("vel = %g %g %g   acc = %g %g %g  _dt = %g\n\n", vel[0][jj],vel[1][jj],vel[2][jj], acc[0][jj],acc[1][jj],acc[2][jj], _dt);
				/*
				if (i==1000) {
					printf("cala_orb_sse()::exit(0);\n");
					exit(0);
				}
				*/
				
				
			
			}
			
			_mm_store_ps(vel[0], vel0jj);
			_mm_store_ps(vel[1], vel1jj);
			_mm_store_ps(vel[2], vel2jj);

			_mm_store_ps(pos[0], pos0jj);
			_mm_store_ps(pos[1], pos1jj);
			_mm_store_ps(pos[2], pos2jj);
			
			//_mm_store_ps((float*)failflag, failflagjj);
			
			
			for (int jj=0; jj<nn; jj++) {
				pos[0][jj] += vel[0][jj]*0.5*dt[jj];
				pos[1][jj] += vel[1][jj]*0.5*dt[jj];
				pos[2][jj] += vel[2][jj]*0.5*dt[jj];
			}
			
			
			for (int jj=0; jj<nn; jj++) {
				if (t[jj]>0) {
					REAL m_over_t(p[J[jj]].Mass/t[jj]);
					REAL *_rho = &rho[J[jj]*dg.n];
					for(int i=0; i<dg.n; i++)  
					if (_rho[i]) _rho[i]*=m_over_t;
				}
				//fclose(__fp[jj]);
				//printf("%i\n", norbit[jj]);
			}
			
			
} 
